clear all;

skyRes = 0.3;
earthRes = 0.2;

skyPos = sinusoidalMap(skyRes);
dOmegaSky = skyPos(:,4);
earthPos = sinusoidalMap(earthRes);

f = linspace(0.1, 8192, 20000);
f0 = min(f);
logLambda = 8.0;
m1 = 10.0;
m2 = 10.0;
psi = 0.0;
iota = 0.0;

   %==============================================================
   j = 1; % (Hanford)
   Sh = detnoisepsd('LIGO', f);
   inspDist = calcinspiraldist( m1, m2, f, Sh, f0); % in Mpc
   %The factor of 2 takes care of optimally oriented binaries (iota = 0.0)
   inspDist = inspDist/2;
   Det = LoadDetectorData('H1');
   [Fp Fc] = ComputeAntennaResponse(skyPos(:,2), skyPos(:,1), psi, Det.d);
   r1 = (inspDist/logLambda) * sqrt((1 + (cos(iota))^2)^2 * Fp.^2 + 4*(cos(iota))^2 * Fc.^2);

   %==============================================================
   j = 2; % (Livingston)
   Sh = detnoisepsd('LIGO', f);
   inspDist = calcinspiraldist( m1, m2, f, Sh, f0); % in Mpc
   %The factor of 2 takes care of optimally oriented binaries (iota = 0.0)
   inspDist = inspDist/2;
   clear Det; Det = LoadDetectorData ('L1');
   [Fp Fc] = ComputeAntennaResponse(skyPos(:,2), skyPos(:,1), psi, Det.d);
   r2 =  (inspDist/logLambda) * sqrt((1 + (cos(iota))^2)^2 * Fp.^2 + 4*(cos(iota))^2 * Fc.^2);
            
   %==============================================================
   j = 3; % (Virgo)
   Sh = detnoisepsd('Virgo', f);
   inspDist = calcinspiraldist( m1, m2, f, Sh, f0); % in Mpc
   %The factor of 2 takes care of optimally oriented binaries (iota = 0.0)
   inspDist = inspDist/2;
   clear Det; Det = LoadDetectorData ('V1');
   [Fp Fc] = ComputeAntennaResponse(skyPos(:,2), skyPos(:,1), psi, Det.d);
   r3 =  (inspDist/logLambda) * sqrt((1 + (cos(iota))^2)^2 * Fp.^2 + 4*(cos(iota))^2 * Fc.^2);


   % Various combinations in double coincidence 
   rmax2.HL = min( [r1 r2]' );
   rmax2.HV = min( [r1 r3]' );
   rmax2.LV = min( [r2 r3]' );
   
   rmax2_coinc_HLV  = max( [rmax2.HL; rmax2.HV; rmax2.LV] );
   coinc2_vol.HLV = (1/3) * sum ((rmax2_coinc_HLV).^3 * dOmegaSky);

   % Various combinations in triple coincidence 
   rmax3.HLV = min( [r1 r2 r3]' );

   rmax3_coinc_HLV  = rmax3.HLV;
   coinc3_vol.HLV = (1/3) * sum ((rmax3_coinc_HLV).^3 * dOmegaSky);


   % Now loop through various locations on the Earth to place the IndIGO detector
   indigoDetPsi = 0; %in degrees
   Sh = detnoisepsd('LIGO', f);
   inspDist = calcinspiraldist( m1, m2, f, Sh, f0); % in Mpc
   %The factor of 2 takes care of optimally oriented binaries (iota = 0.0)
   inspDist = inspDist/2;


   for k=1:size(earthPos,1)

       	clear Det;
       	Det.psi = [indigoDetPsi 00 00]*[1; 1/60; 1/3600];
   	Det.lambda = earthPos(k,1)*180/pi; %degrees
   	Det.phi = earthPos(k,2)*180/pi; %degrees

   	theta = (90-Det.phi)*pi/180;
   	phi   = (Det.lambda)*pi/180;
   	psi   = (Det.psi*pi)/180;

   	thetaCap = ...
        	    [ -cos(theta)*cos(phi); -cos(theta)*sin(phi); sin(theta) ];
   	phiCap = [ -sin(phi); cos(phi); 0.0 ];

   	Det.X    =  cos(psi)*thetaCap + sin(psi)*phiCap;
   	Det.Y    = -sin(psi)*thetaCap + cos(psi)*phiCap;
   	Det.Z    =  cross(Det.X,Det.Y);

   	%----- Response matrix.
  	 Det.d = 0.5*(Det.X*Det.X' - Det.Y*Det.Y') ;

   	[Fp Fc] = ComputeAntennaResponse(skyPos(:,2), skyPos(:,1), psi, Det.d);
   	r4 =  (inspDist/logLambda) * sqrt((1 + (cos(iota))^2)^2 * Fp.^2 + 4*(cos(iota))^2 * Fc.^2);

        % Various combinations in double coincidence 
   	rmax2.HI = min( [r1 r4]' );
   	rmax2.LI = min( [r2 r4]' );
   	rmax2.VI = min( [r3 r4]' );
	
	rmax2_coinc_HLVI = max( [rmax2.HL; rmax2.HV; rmax2.HI; rmax2.LV; rmax2.LI; rmax2.VI] );
	coinc2_vol.HLVI(k) = (1/3) * sum ((rmax2_coinc_HLVI).^3 * dOmegaSky);

        % Various combinations in triple coincidence 
   	rmax3.HLI = min( [r1 r2 r4]' ); 
   	rmax3.HVI = min( [r1 r3 r4]' ); 
   	rmax3.LVI = min( [r2 r3 r4]' ); 

	rmax3_coinc_HLVI = max( [rmax3.HLV; rmax3.HLI; rmax3.HVI; rmax3.LVI]);
	coinc3_vol.HLVI(k) = (1/3) * sum ((rmax3_coinc_HLVI).^3 * dOmegaSky);


   end

  % define the map structure:
   detNames = {'LHO','LLO','Virgo','INDIGO'};
   projectionType = 'mollweid';
   mapstruct = defaultm(projectionType);
   mapstruct.angleunits = 'degrees';
   mapstruct.origin = [0 0 0];
   mapstruct.falseeasting = 0;
   mapstruct.falsenorthing = 0;
   mapstruct.scalefactor = 1;
    
  % load the worldmap coastlines:
   load coast;
   [xc, yc] = mfwdtran(mapstruct, lat, long);

   figure(1); clf; 
   xproject( earthPos, coinc3_vol.HLVI / coinc3_vol.HLV, projectionType );
   hold on;
   plot(xc, yc, 'k-', 'LineWidth', 0.2); % coastlines
   % Plot the locations of all the detectors on top of this
   for j=1:3
   	det = LoadDetectorData( detNames{j} );
   	[x, y] = mfwdtran(mapstruct, det.phi, det.lambda);
   	plot( x, y, 'mo', 'MarkerFaceColor', 'm', 'MarkerSize', 8);
   end
   hold off;
   colorbar;


   
    




   

   






   






